clc; clear;close all;
h = [2/16,2/32,2/64,2/128,2/256];
k=2;
[b_fun, db_fun] = pro_basis();
u0 = @(x) (2/3)+(1/3)*sin(pi*x);
u0_d = @(x) (pi/3)*cos(pi*x);
% u0 = @(x) sin(x);
% u0_d = @(x) cos(x);
errors_L2 = zeros(1,5);
errors_max = zeros(1,5);
[basis,dbasis] = gauss_point_value(b_fun,db_fun);
for d = 3 : 3
    X = -1 : h(d) : 1;
    cfl = 0.18;
    choice = 1;
    N = length(X)-1;
    x_half = 0.5*(X(1:end-1) + X(2:end));
    %生成矩阵
    M = zeros(k+1,k+1);
    % u_true = zeros(1,N);
    for m = 1 : k+1
        for n = 1 : k+1
            PmPn = @(x) b_fun{n}(x).*b_fun{m}(x);
            M(m,n) = (h(d)/2)*quadgk(PmPn, -1,1);
        end
    end
    %对初值做L_2投影
    A0_coef = zeros(k+1, N);
    for i = 1 : N
        x_mid = (X(i) + X(i+1))/2;
        b = zeros(k+1,1);
        for l = 1 : k+1
            Pl = @(x) b_fun{l}((2*(x - x_mid)/h(d)));
            uPl = @(x) u0(x).*Pl(x);
            b(l) = quadgk(uPl,X(i),X(i+1));
        end
        A0_coef(:,i) = M\b;
    end
    dt = cfl*h(d);
    t = 0;
    T = 1.5;
    while t < T
        if t + dt >T
            dt = T-t;
        end
        t = t + dt;
        if choice == 1
            A0_coef = RK3_burger(A0_coef,dt,k,N,M,b_fun,db_fun);
            %A0_coef = TVB_limiter(A0_coef,k,h(d),N,b_fun,basis);
            %A0_coef = minimod_limiter(A0_coef,k,N,b_fun,basis);
            A0_coef = MMP_limiter(A0_coef,k,N,b_fun,basis);
        else
            A0_coef = RK2_burger(A0_coef,dt,k,N,M,b_fun,db_fun);
            %A0_coef = minimod_limiter(A0_coef,k,N,b_fun,basis);
            %A0_coef = TVB_limiter(A0_coef,k,h(d),N,b_fun,basis);
            %A0_coef = minimod_tilde(A0_coef,k,N,b_fun);
            %A0_coef = MMP_limiter(A0_coef,k,N,b_fun,basis);
        end
    end
    u_num = zeros(1,N);
    for i = 1 : N
        for l = 1 : k+1
            u_num(i) = u_num(i) + A0_coef(l,i)*b_fun{l}(0);
        end
    end
    %计算真值
    u_true = zeros(1,N);
    for i = 1:N
        tmp = x_half(i);
        g = @(x) x - (T*(u0(x))+x-tmp)/(T*u0_d(x)+1);
        if x_half(i) < 0
            x0 = 0.49;
        else
            x0 = -0.49;
        end
        b = New_Iter(x0,g,1e-16,2000);
        u_true(i) = u0(b);
    end
    errors_L2(d) = sqrt(h(d))*sqrt(sum((u_true-u_num).^2));
    errors_max(d) = max(abs(u_true-u_num));
end
plot(x_half,u_num,'m*',x_half,u_true,'c-');
legend('数值解','真解');
% orders_L2 = log2(errors_L2(1:end-1)./errors_L2(2:end));
% orders_inf = log2(errors_max(1:end-1)./errors_max(2:end));
% fprintf("inf 误差：");
% for i = 1:length(errors_max)
%     fprintf(" %e", errors_max(i));
% end
% fprintf("\n");
% 
% fprintf("inf 阶数：%s\n", num2str(orders_inf));
% fprintf("L2 误差：");
% for i = 1:length(errors_L2)
%     fprintf(" %e", errors_L2(i));
% end
% fprintf("\n");
% 
% fprintf("L2 阶数：%s\n", num2str(orders_L2));
function x = New_Iter(x0,g,tol,N)
x = g(x0);
n=1;
while abs(x-x0)>=tol && n<N
    x0 = x;
    x = g(x0);
    n = n+1;
end
x = g(x);
end

%生成基函数与基函数导数对应高斯积分点的值
function [basis, dbasis] = gauss_point_value(b_fun,db_fun)
points = [0, -sqrt(3/5), sqrt(3/5)];
%weights = [8/9,5/9,5/9];
basis = zeros(3,3);
dbasis = zeros(3,3);
for l = 1 : 3
    for i = 1 : 3
        basis(i,l) = b_fun{l}(points(i));
        dbasis(i,l) = db_fun{l}(points(i));
    end
end
end
function [b_fun, db_fun] = pro_basis()
b_fun = {};
db_fun = {};
b_fun{1} = @(x) 0.*x +1;
db_fun{1} = @(x) 0.*x;
b_fun{2} = @(x) x;
db_fun{2} = @(x) 0.*x+1;
b_fun{3} = @(x) (3*x.^2-1)/2;
db_fun{3} = @(x) 3*x;
end
function A0 = one_step_burger(A0_coef,k,N,M,b_fun,db_fun)
New_A0 = zeros(size(A0_coef));
[basis,dbasis] = gauss_point_value(b_fun,db_fun);
f = @(x) (x.^2)/2;
% for i = 1 : N
%     %计算每个区间上的近似多项式
%     tmp = @(x) 0.*x;
%     for l = 1 : k+1
%         tmp = @(x) tmp(x) + A0_coef(l,i)*b_fun{l}(x);
%     end
%     uh{i} = tmp;
% end
%存储基函数在-1，1上的值
basis_value = zeros(k+1,2);
for l = 1 : k+1
    basis_value(l,1) = b_fun{l}(-1);
    basis_value(l,2) = b_fun{l}(1);
end
for i = 1 : N
    if i == 1
        q = N;
    else
        q = i-1;
    end
    if i == N
        p = 1;
    else
        p = i+1;
    end
    %计算数值通量
    u_now_right = (A0_coef(:,i))'*basis_value(:,2);
    u_now_left = (A0_coef(:,i))'*basis_value(:,1);
    u_next_left = (A0_coef(:,p))'*basis_value(:,1);
    u_last_right = (A0_coef(:,q))'*basis_value(:,2);
    % u_now_right1 = uh{i}(1);
    % u_now_left1 = uh{i}(-1);
    % u_next_left1 = uh{p}(-1);
    % u_last_right1 = uh{q}(1);
    alpha = max(abs(u_next_left),abs(u_now_right));
    flux_plus = 0.5*(f(u_next_left)+f(u_now_right)-alpha*(u_next_left-u_now_right));
    alpha = max(abs(u_now_left),abs(u_last_right));
    flux_minus = 0.5*(f(u_now_left)+f(u_last_right)-alpha*(u_now_left-u_last_right));
     %flux_plus = (1/2)*(f(point_value(1,i+1))+f(point_value(2,i+1))-alpha*(point_value(1,i+1)-point_value(2,i+1)));
     %flux_minus = (1/2)*(f(point_value(1,i))+f(point_value(2,i))-alpha*(point_value(1,i)-point_value(2,i)));
    %tmp = @(x) uh{i}(x).*uh{i}(x);
    b = zeros(k+1,1);
    int_value = gauss_int_burger(A0_coef(:,i),basis,dbasis,k);
    for l = 1 :k+1
        %tmp1 = @(x) (1/2)*tmp(x).*db_fun{l}(x);
        %a = quadgk(tmp1,-1,1);
        b(l) = int_value(l) -flux_plus*b_fun{l}(1) + flux_minus*b_fun{l}(-1);
        %b(l) = gauss_int(tmp1,2) - flux_plus*b_fun{l}(1) + flux_minus*b_fun{l}(-1);
    end
    New_A0(:,i) = M\b;
end
A0 = New_A0;
end

function int_value = gauss_int_burger(A0_coef,basis,dbasis,k)
%先算u A0_coef k+1*1，basis 3*k+1, dbasis 3*k+1
weights = [8/9,5/9,5/9];
u = zeros(3,1);
u(1) = ((basis(1,1:k+1)*A0_coef)^2)/2;
u(2) = ((basis(2,1:k+1)*A0_coef)^2)/2;
u(3) = ((basis(3,1:k+1)*A0_coef)^2)/2;
int_value = zeros(k+1,1);
for l = 1 : k+1
    tmp = [u(1)*dbasis(1,l);u(2)*dbasis(2,l);u(3)*dbasis(3,l)];
    int_value(l) = weights*tmp;
end
end

function A0_coef = RK2_burger(A0_coef,dt,k,N,M,b_fun,db_fun)
New_A0 = one_step_burger(A0_coef,k,N,M,b_fun,db_fun);
New_A0_coef = A0_coef + dt*New_A0;
New_A0 = one_step_burger(New_A0_coef,k,N,M,b_fun,db_fun);
A0_coef = 0.5*A0_coef + 0.5*New_A0_coef + 0.5*dt*New_A0;
end

function A0_coef = RK3_burger(A0_coef,dt,k,N,M,b_fun,db_fun)
New_A0 = one_step_burger(A0_coef,k,N,M,b_fun,db_fun);
A1_coef = A0_coef + dt*New_A0;
A1_coef = 0.75*A0_coef + 0.25*A1_coef + 0.25*dt*one_step_burger(A1_coef,k,N,M,b_fun,db_fun);
A0_coef = (1/3)*A0_coef + (2/3)*A1_coef + (2/3)*dt*one_step_burger(A1_coef,k,N,M,b_fun,db_fun);
end

function A0_coef = minimod_limiter(A0_coef,k,N,b_fun,basis)
New_A0_coef = zeros(size(A0_coef));
uh_onetilde = zeros(1,N);
uh_twotilde = zeros(1,N);
basis_value = zeros(k+1,2);
for l = 1 : k+1
    basis_value(l,1) = b_fun{l}(-1);
    basis_value(l,2) = b_fun{l}(1);
end
if k == 1
    M = zeros(k+1,k+1);
    for l = 1 : k+1
        M(1,l) = b_fun{l}(-1);
        M(2,l) = b_fun{l}(1);
    end
else
    M = zeros(k+1,k+1);
    for l = 1: k+1
        M(1,l) = b_fun{l}(-1);
        M(2,l) = b_fun{l}(1);
        M(3,l) = (1/2)*quadgk(b_fun{l},-1,1);
    end
end
uh_aver = get_average(A0_coef,basis,N);
for i = 1 : N
    u_now_right = (A0_coef(:,i))'*basis_value(:,2);
    u_now_left = (A0_coef(:,i))'*basis_value(:,1);
    uh_onetilde(i) = u_now_right - uh_aver(i);
    uh_twotilde(i) = uh_aver(i) - u_now_left;
end
for i = 1 :N
    if i == 1
        q = N;
    else
        q = i-1;
    end
    if i == N
        p = 1;
    else
        p = i+1;
    end
    uh_onetilde(i) = minimod(uh_onetilde(i),uh_aver(p)-uh_aver(i),uh_aver(i)-uh_aver(q));
    uh_twotilde(i) = minimod(uh_twotilde(i),uh_aver(p)-uh_aver(i),uh_aver(i)-uh_aver(q));
    u_right = uh_aver(i) + uh_onetilde(i);
    u_left = uh_aver(i) - uh_twotilde(i);
    if k == 1
        b = [u_left;u_right];
        New_A0_coef(:,i) = M\b;
    elseif k==2
        b = [u_left;u_right;uh_aver(i)];
        New_A0_coef(:,i) = M\b;
    end
end
A0_coef = New_A0_coef;
end

function mod = minimod(a,b,c)
if sign(a) == sign(b) && sign(b) == sign(c)
    mod = sign(a)*min(abs(a),min(abs(b),abs(c)));
else
    mod = 0;
end
end

function mod = minimod_tilde(a,b,c,h)
M = (pi^2)/3;
if abs(a) <= M*(h^2)
    mod = a;
else
    mod = minimod(a,b,c);
end
end

function A0_coef = TVB_limiter(A0_coef,k,h,N,b_fun, basis)
New_A0_coef = zeros(size(A0_coef));
uh_onetilde = zeros(1,N);
uh_twotilde = zeros(1,N);
basis_value = zeros(k+1,2);
for l = 1 : k+1
    basis_value(l,1) = b_fun{l}(-1);
    basis_value(l,2) = b_fun{l}(1);
end
if k == 1
    M = zeros(k+1,k+1);
    for l = 1 : k+1
        M(1,l) = b_fun{l}(-1);
        M(2,l) = b_fun{l}(1);
    end
else
    M = zeros(k+1,k+1);
    for l = 1: k+1
        M(1,l) = b_fun{l}(-1);
        M(2,l) = b_fun{l}(1);
        M(3,l) = (1/2)*quadgk(b_fun{l},-1,1);
    end
end
uh_aver = get_average(A0_coef,basis,N);
for i = 1 : N
    u_now_right = (A0_coef(:,i))'*basis_value(:,2);
    u_now_left = (A0_coef(:,i))'*basis_value(:,1);
    uh_onetilde(i) = u_now_right - uh_aver(i);
    uh_twotilde(i) = uh_aver(i) - u_now_left;
end
for i = 1 :N
    if i == 1
        q = N;
    else
        q = i-1;
    end
    if i == N
        p = 1;
    else
        p = i+1;
    end
    uh_onetilde(i) = minimod_tilde(uh_onetilde(i),uh_aver(p)-uh_aver(i),uh_aver(i)-uh_aver(q),h);
    uh_twotilde(i) = minimod_tilde(uh_twotilde(i),uh_aver(p)-uh_aver(i),uh_aver(i)-uh_aver(q),h);
    u_right = uh_aver(i) + uh_onetilde(i);
    u_left = uh_aver(i) - uh_twotilde(i);
    if k == 1
        b = [u_left;u_right];
        New_A0_coef(:,i) = M\b;
    elseif k==2
        b = [u_left;u_right;uh_aver(i)];
        New_A0_coef(:,i) = M\b;
    end
end
A0_coef = New_A0_coef;
end

function A0_coef = MMP_limiter(A0_coef,k,N,b_fun,basis)
New_A0_coef = zeros(size(A0_coef));
max_value = 1;
min_value = -1;
uh_aver = get_average(A0_coef,basis,N);
basis_value = zeros(k+1,2);
for l = 1 : k+1
    basis_value(l,1) = b_fun{l}(-1);
    basis_value(l,2) = b_fun{l}(1);
end
if k == 1
    M = zeros(k+1,k+1);
    for l = 1 : k+1
        M(1,l) = b_fun{l}(-1);
        M(2,l) = b_fun{l}(1);
    end
else
    M = zeros(k+1,k+1);
    for l = 1: k+1
        M(1,l) = b_fun{l}(-1);
        M(2,l) = b_fun{l}(1);
        M(3,l) = (1/2)*quadgk(b_fun{l},-1,1);
    end
end
for i = 1 : N
    %求每个区间上的最大值和最小值
    % x_min = fminbnd(uh{i},-1,1);
    % x_max = fminbnd(@(x)(-1)*uh{i}(x),-1,1);
    % max_inter = uh{i}(x_max);
    % min_inter = uh{i}(x_min);
    [max_inter,min_inter] = get_extrema(A0_coef(:,i),b_fun,k);
    a = (max_value-uh_aver(i))/(max_inter-uh_aver(i));
    b = (uh_aver(i)-min_value)/(uh_aver(i)-min_inter);
    theta = min(1,min(a,b));
    % uh{i} = @(x) theta*uh{i}(x) + (1-theta)*uh_aver(i);
    %算新系数
    u_ori_right = A0_coef(:,i)'*basis_value(:,2);
    u_ori_left = A0_coef(:,i)'*basis_value(:,1);
    u_right = theta*u_ori_right + (1-theta)*uh_aver(i);
    u_left = theta*u_ori_left + (1-theta)*uh_aver(i);
    % u_right = uh{i}(1);
    % u_left = uh{i}(-1);
    if k == 1
        b = [u_left;u_right];
        New_A0_coef(:,i) = M\b;
    elseif k==2
        b = [u_left;u_right;uh_aver(i)];
        New_A0_coef(:,i) = M\b;
    end
end
A0_coef = New_A0_coef;
end

function u_aver = get_average(A0_coef, basis,N)
weights = [8/9,5/9,5/9];
[row,col] = size(A0_coef);
k = row ;
u_aver = zeros(1,N);
for i = 1 : N
    value = basis(:,1:k)*A0_coef(:,i);
    u_aver(i) = (1/2)*(weights*value);
end
end

function [max_value,min_value] = get_extrema(A0_coef,b_fun,k)
if k == 1
    basis_value_right = [b_fun{1}(1),b_fun{2}(1)];
    basis_value_left = [b_fun{1}(-1),b_fun{2}(-1)];
    basis_value = [basis_value_right;basis_value_left];
    value = basis_value*A0_coef;
elseif k==2
    basis_value_right = [b_fun{1}(1),b_fun{2}(1),b_fun{3}(1)];
    basis_value_left = [b_fun{1}(-1),b_fun{2}(-1),b_fun{3}(-1)];
    extrema = (-A0_coef(2))/(3*A0_coef(3));
    if extrema >= -1 && extrema <= 1
        basis_value_extrema = [b_fun{1}(extrema),b_fun{2}(extrema),b_fun{3}(extrema)];
        basis_value = [basis_value_right;basis_value_left;basis_value_extrema];
        value = basis_value*A0_coef;
    else
        basis_value = [basis_value_right;basis_value_left];
        value = basis_value*A0_coef;
    end
end
max_value = max(value);
min_value = min(value);
end